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Abstract 

A method Is presented for computing sensitiv- 
ity derivatives with respect to Independent (Input) 
variables for complex, internally coupled systems, 
while avoiding the cost and inaccuracy of finite 
differencing performed on the entire system anal- 
ysis. The method entails two alternative algo- 
rithms: the first is based on the classical im- 

plicit function theorem formulated on residuals of 
governing equations, and the second develops the 
system sensitivity equations in a new form using 
the partial (local) sensitivity derivatives of the 
output with respect to the input for each part of 
the system. A few application examples are pre- 
sented to illustrate the discussion. The method 
has a potential to answer the "what if" questions 
by presenting engineers with sensitivity informa- 
tion on design trade-offs to guide human judgment 
and formal optimization. In addition, the method 
is compatible with the modern technology of dis- 
tributed computing as well as traditional division 
of design tasks among groups of specialists in the 
design process. The capability to quantify the 
effects of proposed design changes may provide the 
basis for a mathematical model of design. 

Nomenclature 

A.j the i-th CA 

CA contributing analysis, a "black box" 

transforming input into output data used in 
analysis of a system; usually associated with 
an engineering discipline, or a physical part 
of the system 

c Q CPU time for computing system sensitivity 
derivatives by a one-step finite difference 
procedure involving repeated analysis of the 
entire system 

c‘ CPU time for computing system sensitivity 

derivatives using GSE1 

c' ' CPU time for computing system sensitivity 
derivatives using GSE2 

F vector of functions forming the equations 
governing a physical phenomenon 
f functional relationship 

GSE1 Global Sensitivity Equations based on partial 

derivatives of residuals 
GSE2 Global Sensitivity Equations based on the 

partial derivatives of output with respect to 
input of each CA 

H number of input items received by a CA from 

other CA's 

I identity matrix 

M number of independent variables in a CA 

m number of unknown variables in a CA 

N number of the CA's in a system analysis. 

System Analysis one solution of Eq. 1 for all the 
output unknowns Y 

X vector of independent variables 

Y vector of dependent variables 

Z number of unknown variables in a CA 


♦Deputy Head Interdi scipl i nary Research Office, 
AF AIAA. 


0. 0 ,y identifiers for CA's in a small system of 

three CA's, equivalent of A^, A 3 

Subscripts, superscripts, special markings: 

1, j,k subscripts identifying CA's, elements of 

vectors, and elements of matrices 
o superscript or subscript identifying an ini- 
tial value, or a normalization denominator 
overbar linearized function 
tilde normalized, nondimensional quantity 
arrow above character designates a vector 

Other symbols used locally are identified where 
introduced. 

Introduction 

"What if" is the all important question that 
arises again and again in design. Indeed, it may 
be argued that the design process is not complete 
until all such pertinent questions have been asked, 
satisfactorily answered, and the answers translated 
into design changes toward a product as good as it 
can be made under a set of given restrictions. If 
the object being designed is a complex, coupled 
system, the "what if" questions are difficult to 
answer because, to borrow a phrase from (Ref. 1), 
"if you make any change to it there are likely to 
be many subtle consequences". A recent example 
from aerospace is the forebody shape in hypersonic 
aircraft whose change influences structures, aero- 
dynamics, propulsion, and, ultimately, the 
performance. 

Many "what if" questions cannot be quantified 
and engineering judgment is indispensable to answer 
them. However, in aerospace vehicle design a great 
deal of "what if" questions can be quantified ei- 
ther by assessing the effects of relatively large 
variations of the variables involved (a parametric 
study) or by considering very small, theoretically 
infinitesimal variations to calculate sensitivity 
derivatives. 

The focus of this paper is on sensitivity 
analysis. While recent developments in numerical 
methods provided engineers with many useful tech- 
niques for disciplinary, or subsystem, sensitivity 
analysis, e.g., (Refs. 2, 3, 4), examination of 
literature, e.g., (Refs. 1, 5, 6) shows a void as 
far as the comparable methods applicable to entire 
systems are concerned. This paper's purpose is to 
address that void and to offer a system sensitivity 
analysis capable of answering the quantitative 
"what if" design questions. To that end, the paper- 
presents a method for computing sensitivity deriva- 
tives with respect to independent (input) variables 
for complex, internally coupled systems, while 
avoiding the cost and inaccuracy of finite differ- 
encing performed on the entire system analysis. 

The method entails two alternative algorithms: the 

first is based on the classical implicit function 
theorem formulated on residuals of governing equa- 
tions, and the second develops the system sensitiv- 
ity equations in a new form using the partial 
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(local) sensitivity derivatives of the output with 
respect to the input for each part of the system. 

A few application examples are presented to illus- 
trate the discussion. 

Statement of the Problem 

In this paper, a complex, internally coupled 
system is defined as physical object whose behavior 
is described by a vector Y obtainable as a solution 
of a set of simultaneous (coupled) equations which 
can be partitioned into subsets such as 

a((X,Y 0 ,Y Y ),Y a ) = 0 

S((X,Y,Y ),Y ) = 0 (1) 

ay p 

r((X,Y ,Y ),Y ) = 0 
pa y 

Each of the system subsets represents a distinct, 
separate analysis that will be referred to as con- 
tributing analysis (CA), usually associated with a 
particular engineering discipline, or a distinct 
physical part (a subsystem) of the system, or both. 
Partitioning of the system analysis into separate 
but coupled CA's amounts to a system decomposition. 
The Operations Research literature calls such par- 
titioning an aspect decomposition if the CA's cor- 
respond to disciplines, and an object decomposition 
if they correspond to physical subsystems (Ref. 7). 
In most engineering problems both types of decompo- 
sition are used simultaneously to break the large 
task into smaller ones. Mathematics developed in 
this paper applies equally to both types. All the 
mathematical discussion herein is based on three 
partitions because that is a number which is con- 
veniently small, and yet sufficient to establish 
patterns that can easily be generalized to arbi- 
trarily large number of partitions. Solving the 
entire set of equations will be referred to as 
the system analysis which can be written as 
F(Y,X) = 0. 

Each CA yields a solution in form of a vector 
Y* (where subscript ★ stands for a, s or y, and 
identifies a subset of Y) listed last in the paren- 
theses, given the input listed in the inner paren- 
theses. The system is internally coupled because 
the input to one CA includes outputs from the other 
CAs - as shown by the arrows in Figs. 1. The cou- 
pled system is depicted in Fig. 2 by a directed 
graph representation (e.g., (Ref. 1)). 

This paper's focus is on large scale applica- 
tions in which at least some CA's are nonlinear and 
complex, so that the system analysis can only be 
done iteratively. Typically, a CA is carried out 
by a group of specialists, maybe at a separate sub- 
contractor organization. This may be illustrated 
by an example of aircraft wing design incorporating 
nonlinear aerodynamics, structures, and active con- 
trol (aspect decomposition), or substructuring (ob- 
ject decomposition). 

The system solution Y is sensitive to the 
independent variables X present in the CA inputs. 

It is important to emphasize that the independent 
variables X may include not only the designer- 
decided inputs (design variables) but also other- 
inputs external to the system, for example, loads. 


heat flux, etc. In the most general case, all 
variables X may occur in the input to each CA, but 
in most practical applications only a subset of the 
vector X will enter the input of a particular CA. 

One way to compute sensitivity derivatives of 
the solution Y with respect to the independent 
variables X is a finite difference technique de- 
picted by a flowchart in Fig. 3 in its simplest, 
one-step- forward, version. It requires repetition 
of the system analysis for every perturbed X. This 
may be prohibitively costly, particularly if the 
system analysis is nonlinear and/or iterative. 

Even more importantly, it may be inaccurate to the 
point of producing roeanTn^l ess results as the 
effect of small perturbations in X may drown in the 
noise of the iterative solution of the system 
(e.g., Ref. Attempting^ remedy this effect 
by increasing the perturbation magnitude may intro- 
duce significant error due to the analysis nonlin- 
earity. Consequently, the perturbation range in 
which accuracy of finite differencing is acceptable 
becomes problem dependent and may not even exist. 

Thus, the problem is how to calculate the 
sensitivity derivatives of the system solution Y 
with respect to the independent variables X without 
resorting to a finite difference operation involv- 
ing the entire system analysis as in Fig. 3. 

Solution 

As mentioned in the Introduction, there are at 
least two ways of solving the system sensitivity 
problem. A residual -based solution will be intro- 
duced first, and an alternative using local output 
sensitivity will follow. 

Residual -Based Solution 

The implicit function theorem of functional 
analysis, e.g.. Ref. 9 states that a set of 
governing equations 

F ( Y , X ) = 0; Y = f ( X ) (2) 

has the following sensitivity equations 

[s] ($0 ■ • {s,} ,3 > 

The sensitivity equations are always simultaneous, 
linear, and algebraic, regardless of the mathemati- 
cal nature (nonlinear, transcendental, etc.) of the 
governing equations of the system. In Eq. 3, the 
matrix of coefficients, m x m, is a Jacobian matrix 
of the partial derivatives with respect to depen- 
dent variables, and the right-hand- side vector con- 
tains the partial derivatives with respect to a 
particular independent variable. These partial 
derivatives are evaluated using the X and Y values 
which satisfy Eq. 2. In other words, solution of 
the governing equations, Eq. 2, is a prerequisite 
to forming and solving the sensitivity equations, 
Eq. 3. 

The solution vector of Eq. 3 comprises the 
derivatives of the dependent variables with respect 
to a particular independent variable. It will be 
useful in further discussion to have noted at this 
point that Eq. 3 is based on residuals of Eq. 2, 
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i.e., a perturbation of one element in X alone 
would generate a vector of residuals of F replacing 
zero on the right hand side of the equation. Simi- 
larly, a perturbation of one element in Y alone 
would also generate a residual vector. Conse- 
quently, to maintain the right hand side at zero 
despite the perturbation of X, there must be a 
change of in Y subordinated to the change of X to 
make the residual vectors due to Y and X offset 
each other. Equation 3 merely states that to gen- 
erate compensating residuals the rates of change of 
the residuals with respect to the dependent and 
independent variables must balance each other, tak- 
ing into account the implicit dependence of Y on X. 
In other words, the total derivative with respect 
to X of the residuals of Eq. 2 must vanish. 

The method for computing the terms given in 
Eq. 3 is problem-dependent. Obviously, an analyti- 
cal differentiation is preferred but, if that is 
not possible, a finite difference technique may be 
applied. Since the finite difference technique in 
this application is used to calculate the partial 
derivatives of residuals, it requires only an eval- 
uation of F{Y,X) for perturbations of its arguments 
instead of a solution of F(Y,X) = 0 for each per- 
turbation. Thus, the finite difference operation 
performed on the entire system analysis as in 
Fig. 3 is eliminated. 

When applied to the partitioned system in 
Eq. 1, the sensitivity equations 3 take on this 
form 


structural analysis where the residuals are unequ- 
ilibrated loads). These reasons motivated deriva- 
tion of a new form for the system sensitivity equa- 
ions not predicated on the residuals. 

Formulation Based on Sensitivities of Individual 

ZTs 


Residual -independent sensitivity equations may 
be derived in more than one way. The derivation 
shown below is based on linearization of the 
governing equations 1, an alternative derivation is 
shown in Appendix A. 

Equation 1 relate each partition of Y to the X 
and the other partitions of Y so that from each 
equation: 


■ V«'VV 

■ y*-vv 


V - f (X,Y 0 ,Y ) 

Y Y B a 


(5) 


These functions may be linearized in the neighbor- 
hood of the solution of Eq. 1 denoted Y a0 , Y g0 , Y 
using a curtailed Taylor series. Using Y a as an Y 
example: 


= Y n + (3 f „/3 X )AX + (af /3Y UY fl 
a ao a a 0 0 
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»Y.,/>«. |T 
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3 V 3X *k 



+ <3f a /3Yy) A Y Y ; (6) 

By moving all terms to the left hand side, Eq. 6 
are transformed into a linearized version of Eq. 1 


« - Y « - Y a 0 - (3 V 3X)(X ' V 
- ‘ 3 V 3 V (Y B * 'BO 1 

+ - < 3 V 3 V ( \ - v = 0 


« - Y e - V - (9 V 3X >< X - V 


referred to as the Global Sensitivity Equations 1 
(GSE1) . These equations contain as unknowns the 
sensitivity derivatives of the system solution Y 
(partitioned) with respect to an independent vari- 
able X (one at a time). Their matrix of coeffi- 
cients is populated by the partial derivatives of 
the residuals of each CA with respect to the input 
that CA receives from the other CA's, and the 
right-hand-side vector is formed from the partial 
derivatives of the CA residuals with respect to the 
independent variable directly affecting that CA. 

For a general case of N CA*s, the equations acquire 
a format shown in Appendix A. 

Despite their potential cost and accuracy 
advantages, the use of the sensitivity equations 4 
based on residuals may not be straightforward in 
engineering practice because existing disciplinary 
codes have usually no provisions to compute the 
residuals, and the residuals usually have no obvi- 
ous physical meaning that would allow the user to 
judge validity of the numbers (An exception is 


• - '.o' 

* - (lf„/)Y ) (Y - Y ) ■ 0 

0 Y Y YO 


y - \ - \ 0 - (3 V 3X)(X - V 

- (3f /aY ) (Y - Y ) 

Y a a aO 


+ - < 9 V 3 V< y b - V = 0 


(7) 


Under some conditions discussed in Appendix A, 

Eq. 7 may have no solution because of singularity 
of their matrix of coefficients. Assuming that 
singularity conditions examined in the Appendix A 
do not occur, it is axiomatic that Eq. 7 and 1 have 
the same solutions Y and that these solutions have 
the same derivatives with respect to X, consequent- 
ly, we may treat Eq. 7 as surrogate governing equa- 
tions. The implicit function theorem may now be 
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applied to these equations just as it was applied 
to Eq. 1 by performing the differentiation shown in 
Eq. 3. This yields sensitivity equations in the 
form: 


I 

' -3f /3Y ’ 

1 a 3 1 

- 3 V 3 \ 

‘ 3 V 3Y a 

p a 

1 1 
1 I 

-»v*\ 

-3f /3Y 
Y <* 

: " 3 V 3v e : 

I 

r 'y 

r -* 



3Y a /3X k 


K /a v = s 3f * /3 v 


L 3 V 3X kJ 


l °, 


( 8 ) 


termed Global Sensitivity Equations 2. For a 
general case of N CA's, the equations acquire a 
format shown in Appendix A. 

Equation 8, contains no residuals of the CA's. 
Instead, its matrix of coefficients is populated by 
the sensitivity derivatives of each CA output with 
respect to that CA's input, and its right hand side 
vector represents sensitivity of a CA's output with 
respect to the independent variable (one at a time) 
directly affecting that CA. As far as the complete 
system is concerned, these derivatives are partial 
(local) derivatives, while the solution of Eq. 8 
yields the derivatives of the solution Y 
(partitioned) with respect to an independent 
variable (one at a time). By definition, the 
partial derivatives represent sensitivity of each 
isolated CA, and the derivatives of Y represent the 
system sensitivity with all the couplings (e.g., 
Figs. 1 and 2) fully accounted for. 


Both GSE1 and GSE2, Eqs. 4 and 8, produce the 
same solution vector and both are equally exact 
because they are derived from a mathematical theo- 
rem without any simplifying assumptions or approxi- 
mations. Similarity of mathematical characteris- 
tics and potential usage among the two sets of 
equations allows to limit the ensuing discussion to 
GSE2, occasionally using the notation GSEx to 
address both GSE1 and GSE2. This emphasis on GSE2 
does not imply an unqualified recommendation of 
GSE2 over GSE1. The choice is up to the user and 
it depends on the factors already stated in the 
foregoing as motivations for the GSE2 development, 
and on the considerations of cost and benefits 
discussed later (supported by Appendix B). 


Figure 4 illustrates details of the GSE2 
structure for three CA's. The general pattern is 
easily extrapolated to larger number of CA's. The 
matrix of coefficients has identity submatrices on 
the diagonal. Each off-diagonal submatrix is a 
Jacobian matrix corresponding to a CA. For ex- 
ample, the Jacobian at the upper right corner in 
Fig. 4 contains partial sensitivity derivatives of 
every item of output from the a CA, a column m 
long, with respect to every one of m^ input items 
the a CA receives from the y CA, hence the m^ 


columns in the Jacobian. The corresponding 
Jacobian at the lower left corner comprises the 
partial derivatives of the output from the y CA 
with respect to the input the y CA receives from 
the a CA. In general case the two Jacobi ans are 
not symmetric. An example of the above two 
Jacobian matrices might be drawn from a case of an 
actively controlled flexible wing. Then, assuming 
the CA's a, 6, and y to be aerodynamics, struc- 
tures, and active controls, respectively, the upper- 
right Jacobian would include the partial sensitiv- 
ity derivatives the aerodynamic pressure at select- 
ed locations on the wing to the control surface 
deflections. Correspondingly, the lower left 
Jacobian might comprise of the partial derivatives 
of the control surface deflections to the aerodyna- 
mic pressure coefficients. These derivatives de- 
rive from the control law that establishes a func- 
tional relationship between the control surface 
deflections and the aerodynamic pressure on the 
wing (sensed directly or indirectly). 

On the right hand side, there is a vector of 
the partial sensitivity derivatives of the CA 
output with respect to a particular independent 
variable X. These partial derivatives are nonzero 
for those CA's that are directly influenced by that 
particular independent variable. Referring again 
to the above example of a flexible wing, if the 
variable X were, say, the planform aspect ratio, 
the nonzero elements of the right-hand-side vector 
would occur at the locations corresponding to the 3 
and y partitions, since only the aerodynamics and 
structures would be directly influenced. 

The GSE2 matrix of coefficients depends only 
on the coupling among the CA's and not at all on 
the sensitivity to the design variables. The 
opposite is true for the right-hand- side vectors. 
Thus, the matrix can be formed and factored once, 
and the solutions for many design variables can be 
obtained by repeated back-substitutions of each 
right-hand-side vector. The coupling among the 
CA's is reflected in the topology of the GSE2 
matrix as shown in a few examples in Fig. 5 and 6. 
If there are no couplings (#1), the matrix is an 
identity matrix and the derivatives of Y are equal 
directly to the partial derivatives on the right- 
hand-side. Each coupling link generates an off- 
diagonal Jacobian until the matrix becomes fully 
populated for a fully coupled system (#7). The 
coefficient matrix in the GSE2 exhibits a the same 
pattern. 

There is a coincidence of form between the 
matrix of coefficients in GSEx and the so-called 
equation precedence matrix (or N-Square Matrix) 
used in Operations Research literature (e.g., 

(Ref. 1, 0 . 87)) to analyze internal couplings in 
systems. Namely, each non-zero, off-diagonal 
Jacobian in the GSEx matrix of coefficients 
corresponds to a non-zero element in the N-Square 
matrix for the same system. 

As a matter of a particular interest to a 
structural engineer, one may observe that if F in 
Eq. 2 represents structural load-deflection equa- 
tions, then the Jacobi ans in Eq. 8 correspond to 
substructures or, ultimately, individual finite 
elements. 

In some applications it may be convenient (for 
instance, when the X variables are measured in 
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different units) to have all tertns in GSE2 dimen- 
sionless. A nondimensional version of the GSE2 is 
given in Appendix A. 

Exampl es 

Since both the GSE1 and GSE2 (Eqs. 4 and 8, 
respectively) are rigorously derived from a fun- 
damental theorem they do not need numerical verifi- 
cation. However, a few examples are provided for 
the usage of GSE2 to support the discussion of 
costs and benefits that will follow. 

A simple example of a 2D airfoil in airflow is 
shown in Fig. 7. The airfoil is supported by two 
linear springs attached to a ramp whose angle of 
inclination $ is an independent variable. The 
elastic degrees of freedom allowed are only the 
pitch and the plunge. The lift coefficient is 
assumed to be a nonlinear function of the angle 
of attack illustrated in Fig. 8 and defined in 
Table 1 - Aerodynamics. The function is set up 
deliberately as a transcendental function to admit 
only an iterative system analysis. The angle of 
attack e depends on the ramp angle (design variable 
if, and the airfoil elastic support pitch angle <f>. 

The airfoil on springs is an aerodynamic- 
structure system abstracted as a directed graph in 
Fig. 9. All the equations that constitute the 
Aerodynamic and Structures CA's in the graph are 
listed in Table 1 which also shows the problem 
notation and its correspondence to the generic 
notation used in the paper, and the numerical data 
for the example. The purpose of the example is to 
show computation of the derivatives of the system 
solution output - the lift L and the elastic pitch 
angle $ - by means of the GSE2 and to compare the 
results with those from a finite difference 
technique. 

The system solution was found iteratively and 
is listed in Table 2 for arbitrary ^ value of .05 
rad. Next, the sensitivity derivatives of L and ^ 
with respect to the angle ^ were obtained by the 
finite differences procedure at the system level 
illustrated in Fig. 3 which required repetition of 
the iterative solution for the angle if, incremented 
by .0025 rad to .0525 rad. These derivatives are 
shown in Table 2 and provide reference for compari- 
son with the same derivatives computed using the 
GSE2. Incidentally, the derivative of L is greater 
than the partial derivative due to the elastic 
effect. The GSE2 and the numerical values of the 
partial derivatives that enter these equations are 
also given in Table 2 (these equations are also 
shown in a dimensionless format in Appendix A). 

The partial derivatives were obtained by the same, 
simple, one-step-forward, finite difference pro- 
cedure referred to above but appl ied separately to 
Aerodynamics and Structures CA^s. Fi nally , Table 2 
presents the 6SE? solution that agrees with the 
finite difference results obtained at the system 
level . 

The second example shows how the GSE2 equa- 
tions for a system are made up of the partial 
derivatives for the system CA's. The system is a 
flexible wing with an active control intended to 
reduce the root bending moment. The system di- 
rected graph and the coupling information are shown 
in Fig. 10 - upper part. The bottom part of the 
figure illustrates the make-up of the GSE2. 


Dimensions of the arrays entering the GSE2 
depend on the number of the individual pieces of 
data (coupling channel bandwidth, referred to 
as bandwidth, for short) communicated from one CA 
to another. These dimensions have a strong impact 
on the computational cost of the method as shown in 
the next section and in Appendix B, therefore, it 
is important to keep the bandwidths as small as 
possible. In this example, the Structures-Acti ve 
Control channel does not need to transmit more than 
a few strain gage readings. Similarly, the 
Aerodynamics-Active Control channel transmits only 
a few dynamic pressure sensor indications (or only 
the Mach number and the angle of attack value from 
which the pressures may be inferred) and one, or 
two, control surface deflection angles. In con- 
trast, the information moving along the 
Aerodynamic-Structures channel may include hundreds 
of the dynamic pressure values for discrete loca- 
tions on the wing, if a panel-based CFD code is 
used, and thousands of the nodal point displace- 
ments output from a finite element code. It is 
evident, that this channel will require attention 
to reduce its bandwidth. Such reduction rnay be 
achieved by representing deformations and loads by 
a relatively small number of generalized coordi- 
nates and corresponding generalized forces based on 
modal analysis, following the practice well estab- 
lished in aeroel astici ty analysis. 

Another example of the use of the GSE2 for a 
system with active control is described in 
(Ref. 8). 

Costs and Benefits 

By using GSEx (Eqs. 4 or 8), the cost of 
repetitive system analysis required by finite 
difference procedure (Fig. 3) is eliminated, but 
the cost of generating the input into these equa- 
tions and solving them is added. Using the CPU 
time as a simplified measure of the computational 
cost, Appendix B shows, under a set of assumptions, 
that the cost for the finite difference procedure 
of Fig. 3 increases with the square of the number 
of CA's in the system. On the other hand, the cost 
of generating the input into GSEx under the same 
set of assumptions is proportional to the product 
of the number of CA's and the bandwidth. These 
relations suggest that there is a limit on the 
finite difference procedure's applicability to 
large systems, and show the importance of the band- 
width to the cost of the GSEx. As far as the 
accuracy of GSEx is concerned, it depends on the 
conditioning of its matrix of coefficients (see 
Appendix A) but is not affected by the system 
dimensional ity. 

The principal qualitative advantage of the 
sensitivity analysis based on the GSE2 is that it 
allows to treat the system as decomposed into a set 
of "black boxes" coupled by a well-defined sets of 
data. Each black box may, then, be subjected to 
its own sensitivity analysis performed by special- 
ists intimately familiar with the specifics. The 
specialists may use any means for the partial sen- 
sitivity analysis available such as: finite dif- 

ference procedures, historical statistical data, 
approximate methods, or even judgmental assessment. 
It should also be stressed that they may also draw 
on the disci pi inary , quasi-analytical sensitivity 
analysis algorithms that are now undergoing an in- 
tensive development (Ref. 6). They may even obtain 
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the sensitivity data experimental ly . In general, 
the approach divides the labor and thus creates 
opportunity for concurrent data processing in the 
contemporary distributed computing environment, and 
supports a broad work front in the engineering 
organization. 

Another benefit from the "black box" approach 
is that the GSEx are inherently recursive, in the 
sense that each of the system's "black boxes" (the 
CA's) may be a complex system within itself. If 
so, its sensitivity analysis may be carried out as 
described herein, to produce the sensitivity deriv- 
atives that will be treated as the sensitivity par- 
tial derivatives in the GSEx of the parent system. 

Regarding of the choice of GSE2 vs. GSE1, if 
the computational cost was the only factor (see 
Appendix B) GSE1 would be recommended over GSE2. 
However, the nonavailability of the residuals in 
existing disciplinary codes, and difficulties with 
physical interpretation of the residuals clearly 
favor the GSE2 format. Furthermore, the disciplin- 
ary sensitivity analyses are formulated to yield 
data compatible with input to GSE2 but not GSEI. 
These considerations may be overridden in the fu- 
ture by new disciplinary code developments (avail- 
ability of the residual options) encouraged by the 
strong cost advantage of the GSEI. For now, the 
choice is judgmental. 

Usage in Design 

Systematic procedure for generating the system 
sensitivity data in a design process using the GSE2 
may be organized in a way shown by a Chapin-format 
(Ref. 10) flowchart in Fig. 11. It begins with the 
system analysis for a given X. The partial sensi- 
tivity derivative computations follow in each CA 
independently for the given X and given Y (the 
latter obtained from the system analysis). The 
partial derivative calculations may be carried out 
concurrently . The final operation is an assembly 
and solution of the GSE2. The usage of GSEI is 
similar. The results of a system sensitivity anal- 
ysis may be used to identify the "design drivers", 
and to determine design modifications toward im- 
provement, or they can be input into a formal opti- 
mization procedure. Since, in general case, the 
system solution and its sensitivity analysis have 
to be updated after moving away from the previously 
solved design, the flowchart procedure has to be 
iterated as the design process advances. 

The GSEx may also be used to assess the coup- 
ling strength between any two CA's. This can be 
done by computing derivatives of the GSEx solution 
to the elements of the GSEx matrix in a manner- 
shown below for the GSE2 used as an example. 

The GSE2 may be written as 

[A] {^k} = {RHS}: (9> 

where A and RHS are the matrix of coefficients and 
the right hand side vector, respectively, defined 
in Eq. 8. Since the equations are linear, the 
derivatives of their solution with respect to 
the elements of the matrix A may be obtained by 
substituting Eq. 9 for F in Eq. 2 and, then, 
writing the corresponding sensitivity equations 
using the differentiation pattern of Eq. 3 



In this set of equations the matrix aA/aA^j is all 
empty except unity at the location occupied by the 
element Ay in the matrix A. The vector of the 
partial derivatives of Y with respect to X k is 
available from the solution of Eq. 9, and the deri- 
vatives of the RHS in Eq. 9 with respect to Ay are 
null so they do not appear in Eq. 10. Consequent- 
ly, the unknown derivatives a (a Y/a X k ) /aA i j may be 

obtained by backsubstitution of the new right hand 
side vector over the matrix A decomposed once in 
the solution of Eq. 9 and saved. 

These derivatives measure the influence of the 
partial sensitivity derivatives Ay on the sensi- 
tivity of the system with respect to X and may be 
adopted as indicators of the strength of the coup- 
lings among the parts of the system. A full survey 
of the coupling strengths would require solution of 
Eq. 9 and 10 for each combination of Ay and X^. 

In case of Y and X expressed in nonhomogeneous 
physical units, a dimensionless form shown in 
Appendix A for Eq. 8 would have to be used to ob- 
tain the coupling strength indicators that could be 
compared with each other. The comparison would be 
useful to identify relatively weak couplings that 
might be dropped from the system's mathematical 
model. Thus, the coupl i ng strength indicators may 
augment the system analyst's judgment in searching 
for a compromise between the system model simplic- 
ity and its predictive accuracy. 

Conclusions 

The paper addresses problem of a sensitivity 
of a complex, internally coupled system behavior 
(response) to changes in independent variables. It 
is assumed that the system analysis is made up of 
self-contained analyses, correspondi ng to disci- 
plines and/or physical subsystems, which exchange 
input/output data with each other. The problem is 
solved by formulating rigorous sensitivity equa- 
tions - the global sensitivity equations - derived 
in a form based on the residuals of the governing 
equations, and in a new form that does not rely on 
the residuals. The latter is judged to be more 
useful to engineers than the former. It allows 
evaluation of system sensitivity to independent 
variables on the basis of the partial derivative 
information obtained locally within each contri- 
buting engineering discipline, or within each 
physical subsystem analysis, consistent with the 
decomposition of the design process among the 
specialty groups, and compatible with the tech- 
nology of distributed computing. The equations 
eliminate the need for costly and potentially 
inaccurate finite differencing performed on the 
entire system analysis, and are capable of accept- 
ing experimentally obtained sensitivity data. 

Their computational cost advantage over the 
reference finite difference procedure increases 
with the number of self-contained analyses into 
which the system analysis can be partitioned, and 
is reduced proportionally to the volume of the 
coupling information. 

Derivatives of the solution of the sensitivity 
equations with respect to their own coefficients 
may be useful as indicators of the strength of the 
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couplings among the parts of the system. Ranking 
these indicators by their magnitudes may identify 
weak couplings that might be eliminated from the 
system's model to make it simpler without sig- 
nificant loss of its predictive accuracy. 

The global sensitivity equations in either 
form are offered as a tool to support the design 
process by contributing the system sensitivity in- 
formation as an aid for human judgment and/or for 
use in formal optimization. Inasmuch as the global 
sensitivity equations quantitatively answer, with 
the first order of the accuracy, the "what if" 
questions underlying the design process, they may 
be regarded as a first order mathematical model of 
that process. 
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Appendix A 

Alternative Derivation, General ization to a 
Case of N CA's, and a Dimensionless Format for the 
Global Sensitivity Equations 2. 

Alternative Derivation 


An alternative derivation begins with Eq. 5, 
differentiated with respect to one particular- 
independent variable, say, X k . Using the chain 
rule, the derivatives of Y are 

9 V 3X k = (3f a /3X k } + l3f a /3Y B )(3Y 6 /3X k ) 

+ (af /aY v )(3Y v /9X. ) 

3 V 9X k = (3 V 3X k> + (3 V 3Y a )(3 V 3X k> 

+ (3 /3 Y ) (3 Y /a X.) 

P T 7 * 

3Y Y /3X k = (3 f /a X k ) + (3f Y /3Y a )(3Y a /3X k ) 

+ (3f Y /3Y g )(3Y g /3X k ) (Al) 


Collecting the given and unknown terms, and re- 
arranging, yields the GSE2 in the format of Eq. 8. 

Possible Singularity of the Matrix of Coefficients 
in tq. 7 and GSEx . 

When a solution of a set of nonlinear equa- 
tions, corresponding to Eq. 2, exists by virtue of 
intersection as illustrated for an example of two 
functions in a two-variable space in Fig. 12(a), 
then, the same solution coordinates are defined by 
a set of corresponding linearized equations, 
corresponding to Eq. 7, and represented by the 
dashed tangents in the figure. However, in case of 
the solution existing by virtue of tangency, as in 


Fig. 12(b), the corresponding tangents overlap. 
Hence, their equations, Eq. 7, become singular and 
have no solution. Therefore, they no longer can be 
used as a substitute for the original, nonlinear- 
equations, or as a basis for deriving linear sensi- 
tivity equations, since the matrix of coefficients 
in these equations will also be singular. It will 
be so because Eq. 7 may be written as 

CA(X)]{Y(X)1 = { RHS(X)} (A2) 

and the corresponding sensitivity equations ob- 
tained either from Eq. 4 or Eq. 8 are 


[A(X)]{3Y(X)/3X> = -[3A(X)/3X]{Y> 

+ {3 RHS(X ) /3 X) ; (A3) 

Since Eq. A2 and A3 share the same matrix of coef- 
ficients, its singularity affects both sets of 
equations. 

One may add that the tangency-type solution to 
Eq. 2 depicted in Fig. 12b has also a drawback of 
rendering the coordinates of the solution point S 
ill -defined. As shown in Fig. 13, due to error in 
numerical definition of the tangent functions, the 
coordinates of the point of tangency fall into 
broad intervals of uncertainty. This may also 
occur for the intersection type-type solution, if 
the intersection angle is small, then the matrix of 
coefficients A ( X ) may be ill-conditioned, although 
non-singular. Occurrence of these cases in design 
usually indicates that the design analysis is ill- 
posed and should be reformulated. 

The GSEx matrix of coefficients may also be 
singular, if the system is physically unstable, 
e.g., a wing divergence. However, such instability 
would normally manifest itself at the prerequisite 
stage of the system analysis (solution of Eq. 1), 
so it is not expected to become a problem in the 
sensitivity analysis. 

GSEI and GSE2 Generalized to a Case of N CA's. 


A simple extrapolation of the pattern from the 
case of three CA's discussed in conjunction with 
Eq. 1, 4, and 8, leads to Eq. 1, GSEI, and GSE2, 
respectively , taking on the form of Eq. A4, A5, A6: 


A-( (X 


V- Y 1->>V = 0 


i = 1.N 
j * i;l * i 


{M) 


[B] 


[B] 


{3 Y/a X k > = - {3 A i /3 X^} 

3 V 3Y j 

{3V/3X k } = (3f 1 /3 X k > 

- 3 V 3Y r B ii = 1 


(A5) 


(A6) 
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Di mensionless Form 


To transform the GSE2 from the form of Eq. 8 
to a dimensionless form, we normalize the variables 
Y and X in Eq. 1 to unity by dividing them by their 
initial values (if any of them is zero, a suitable 
nonzero value is used instead). The normalization 
yields, showing two partial derivatives as examples 


OY ai /3Y 6j ) = OV cti /3Y Bj )(Y« j /Y° i ) 


(A7) 


■ OV 3 X k )(XO k/Y°i) <A 8 > 


This normalization introduced into Eqs. 5, 6, and 7 
leads to a dimensionless form for Eq. 8. For in- 
stance, the partition $ of Eq. 8 becomes 


3f 6 
' w 

a 


3f « 3Y a 3f s 

■I 1 T ' P 1 ‘ P 1 * P 1 

8ai j* 1 ■" W~ 8xki 37^ 3xk 


( A9 ) 


where: 


3. repeated solution of the system for a small 
perturbation of a design variable requires 
P2 < p^ iterations. 

4. each CA is directly influenced by M 
independent variables so that the total 
number of independent variables is MN. 

5. there are Z unknown variables Y in each CA. 

6. each CA receives H input variables from the 
remainder of the system. 

7. partial sensitivity derivatives of a CA are 
computed by finite differences, and each CA 
solution repeated for a small perturbation 
of input requires C2 < c^ of CPU time. 

8. computation of residuals for a CA requires 
c 4 CPU time. 

9. solution of the GSE2 and GSE1 with multiple 
right hand sides takes C3 and C5 CPU times, 
respectively. These times are expected to 
be relatively small, due to the use of the 
parallel and vector processing technology. 


0aij 


Y° 7Y . • 1 
al 7 BJ* 3Y 1 J 


W 


X k' Y B °J 


Referring to the example defined in Tables 1 and 2, 
the GSE2 shown in Table 2 transform to the follow- 
ing dimensionless form: 


Under this assumption, the cost of the finite 
difference procedure is a sum of two terms: a 

reference system analysis, and the system analysis 
repeated for small perturbation of each of the in- 
dependent variables. The two terms form, respec- 
tively, the following expression: 


Y e L; Y = *; L = 502.2995 M; 
o fci 


Cq - Nc^p^ + Nc jpgMN 


4 = .0175805 rd 


C 1 P 1 ( N + ^ P 2 /p 1 1Mn2 


(Bl) 


1 

1 


= y ; f = .05 rd 


-.343 

1 


(alT/df^ 

H/rfJ 


r-3431 

0 


(A10) 


in which the great numerical disparity of the terms 
visible in the dimensional form of the equations in 
Table 2 is avoided which is one benefit from the 
dimensionless form. 


For the GSE2-based sensitivity analysis the cost is 
a sum of the cost of one system analysis, the cost 
of solution of each of the N CA's repeated for a 
small perturbation of its (H + M) inputs to compute 
the partial derivatives, and the cost of solving 
the GSE2. These three cost contributions are re- 
presented by the respective terms in the following 
expression: 

c" = Nc 1 P 1 + c 2 N (H + M) + C 3 (B2) 


Appendix B 

Computational costs of the GSEx-based sensi- 
tivity analysis and the reference finite difference 
procedure measured by the CPU times are influenced 
by a very large number of the problem-dependent and 
computer type-dependent (hardware and software) 
variables, but a reasonable estimates can be made 
if a number- of simplifying assumptions are intro- 
duced. The assumptions used here are: 

1. there are N CA’s, each having the same CPU 
time C| for one solution. 

2. complete solution of the system is itera- 

tive; it requires pi iterations, and Nc^ CPU 

time in each iteration. 


Finally, for the GSEl-based sensitivity analysis 
the cost is a sum of the cost of one system anal- 
ysis, the cost of computing the residuals of each 
CA for small perturbations of each of: its unknown 

Y variables, each of the coupling Y variables, and 
each of the X variables influencing that CA. The 
total cost includes also the cost of solving the 
GSE1. These three cost contributions are repre- 
sented by the respective terms in the expression: 

c' = N 1 p 1 + c 4 (Z + H + M) + c 5 (B3) 

The above equations reveal that the finite differ- 
ence procedure cost may be tending toward over- 
whelming values for large number of CA's because of 
the presence of the N 2 * * term in Eq. Bl. On the 
other hand, the cost of the GSEx-based sensitivity 
analysis does not depend on N 2 but is proportional 
to H and the cost of solution of the GSEx. This 
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suggests that its cost advantage over the finite 
difference procedure will increase with the size of 
the system measured by N, provided that magnitude 
of the coupling bandwidth H is judiciously kept 
under control, and the full advantage is taken of 
the progress in computing technology to keep c 3 and 
c 5 as low as possible. 

It is apparent from Eq. B2 and B3 that the 
GSE1 cost has a potential of being much smaller 
than that of GSE2 because evaluation of the 
residuals for a CA takes much less time than 
computation of the partial derivatives of its 
output with respect to its input. Therefore, one 
may expect c^ << c 3 - a reasonable estimate would 
be a one, or two, orders of magnitude difference. 
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TABLE 1. DEFINITION OF EXAMPLE 1 
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TABLE 2. SENSITIVITY ANALYSIS OF EXAMPLE 1. 


SYSTEM SOLUTION 


L = 502.3 N; « = .0176 rd 


DERIVATIVES WITH RESPECT TO Y BY FINITE 
DIFFERENCES 


AY = .0025 

= 14925.16 M/rd; || = .5221287 rd/rd 


GSE2 


Symbolically: 

I -af a /34> 

|3L/3y| 


-3f @ /3L/ I 

(3<|>/3yJ 

Numerically: 

1 -9805.105 

(3L/3y| 


-.3510" 4 1_ 

[34/3yJ 



DERIVATIVES WITH RESPECT TO Y FROM GSE2 


14928.12 N/rd; ||= .5224841 rd/rd 

or OT 


^Numerical values of af a /3<t> and 3f a /3t are 
equal because of relation marked # in Table 1 



Yy +AYy . 

Fig. 3 Finite difference procedure involving the 
system analysis. 


I i 

! aY a /aX K| . 

! a v®M ' . 

“ 9 f y j / 3 Y a i \ / \ 

Fig. 4 Anatomy of the GSE2. 




Bf V_ j 

1J 



Fig. 1 Function vectors (a) forming a set of 
coupled equations (b). 



a (3 Y 4 a P Y 



Fig. 5 System couplings reflected in the GSE2 
matrix. 



Fig. 2 Directed graph representation of the 
system shown in Fig. 1. 





Fig. 6 More examples of system couplings and their 
reflection in the GSE2 matrix. 
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System 


a ( ) = 0 

f|TT-=-o 

Iv ( )=0 

| Solve the system 


GSE2 — 


Use the total derivatives to 
redesign toward improvement 


Fig. 7 Example 1: a simple aerodynamic- 

structures system 



1 O 0 

o 

Fig. 8 Nonlinear relationship C L vs. angle of 
attack in Example 1. 


Fig. 11 Sensitivity procedure in design process 


t,=o v* 

"ll 

//A=° 


hy 0 
// f 2 = 0 


^f 2 = 0 


Figure 12.- Constrained minimum defined by 

intersection (a) and tangency (b). 



Fig. 9 The system from Fig. 7 abstracted as a 

"black box" with a directed graph showing 
internal coupling. 


J 2 = 0 


, f ,=0 




Wing sy stem^aerodynamics 4- structures + active controls 


Figure 13.- Ill-defined constrained minimum. 


fiv 


stress due to r root 
bending moment V 


r-Sensitivity of loads to struct, deform. 

/ /—Total sensit. of loads to struct. 

"i (Ml [o' 1 f o ' sizin 9 


p HIM 




-Partial sensitivity of 
control surface 
deflections to control 
law parameters. 


Total sensitivity of / ,aw parameters. 

deformation to / „ . , 

struct, sizing. ^Part.at sensitmty of 

deformation to struct, sizing 


Lift L = / loads dS;dL/d(str. siz.) = / d(loads)/d($tr. siz.)dS: 


Fig. 10 Example of a flexible wing: a system 

comprising aerodynamics, structures, and 
active control. 
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